FigureS3_H2O2

Load Data

dataset <- read.delim("raw_data/SFigure3_H2O2.txt", stringsAsFactors = FALSE)

dataset$genotype <- gsub(" ", "", dataset$genotype )
dataset$genotype <- factor(dataset$genotype)
dataset$Experiment <- factor(rep(paste0("exp", 1:(length(dataset$genotype)/length(levels(dataset$genotype)))), each=length(unique(dataset$genotype))))

dataset$siRNA <-  factor(gsub(".*[T,O]\\+","",dataset$genotype))
dataset$genotype <-  factor(gsub("\\+.*","",dataset$genotype))

dataset$UID <- factor(paste(dataset$Experiment, dataset$genotype, dataset$siRNA))
dataset$GSID <- factor(paste(dataset$genotype, dataset$siRNA))

# wide format
kable(dataset, row.names = F)

genotype	NT	H2O2_5uM	H2O2_30uM	H2O2_50uM	Experiment	siRNA	UID	GSID
WT	2420	1528	977	248	exp1	siCtrl	exp1 WT siCtrl	WT siCtrl
WT	2180	755	303	60	exp1	siPARP1	exp1 WT siPARP1	WT siPARP1
ALC1KO	2428	923	244	95	exp1	siCtrl	exp1 ALC1KO siCtrl	ALC1KO siCtrl
ALC1KO	2188	834	296	74	exp1	siPARP1	exp1 ALC1KO siPARP1	ALC1KO siPARP1
WT	2072	1332	598	194	exp2	siCtrl	exp2 WT siCtrl	WT siCtrl
WT	2504	624	344	111	exp2	siPARP1	exp2 WT siPARP1	WT siPARP1
ALC1KO	3160	1222	388	122	exp2	siCtrl	exp2 ALC1KO siCtrl	ALC1KO siCtrl
ALC1KO	3536	1296	396	33	exp2	siPARP1	exp2 ALC1KO siPARP1	ALC1KO siPARP1
WT	1768	896	562	219	exp3	siCtrl	exp3 WT siCtrl	WT siCtrl
WT	1888	696	422	160	exp3	siPARP1	exp3 WT siPARP1	WT siPARP1
ALC1KO	2455	712	499	140	exp3	siCtrl	exp3 ALC1KO siCtrl	ALC1KO siCtrl
ALC1KO	3222	1224	467	67	exp3	siPARP1	exp3 ALC1KO siPARP1	ALC1KO siPARP1

library(reshape2)
# reshape to long format
dataset <- melt(dataset, variable.name = "Treatment", value.name = "Counts")

dataset$genotype <- relevel(dataset$genotype, ref = "WT")
dataset$siRNA <- relevel(dataset$siRNA, ref = "siCtrl")
dataset$UID <- relevel(dataset$UID, ref = "exp1 WT siCtrl")

dataset$H2O2 <- gsub("NT","1",dataset$Treatment)
dataset$H2O2 <- gsub("H2O2_|uM","",dataset$H2O2)
dataset$H2O2 <- log10(as.integer(dataset$H2O2))




dataset$Offset <- NA
for(uid in levels(dataset$UID)){
        dataset$Offset[dataset$UID == uid] <- mean(dataset$Counts[dataset$UID == uid])
}

dataset$NormCounts <- dataset$Counts / dataset$Offset



dataset$Offset2 <- NA
for(gsid in levels(dataset$GSID)){
        dataset$Offset2[dataset$GSID == gsid] <- mean(dataset$NormCounts[dataset$GSID == gsid & dataset$H2O2 == 0])
}

dataset$NormCounts2 <- dataset$NormCounts / dataset$Offset2



# long format
kable(dataset, row.names = F)

genotype	Experiment	siRNA	UID	GSID	Treatment	Counts	H2O2	Offset	NormCounts	Offset2	NormCounts2
WT	exp1	siCtrl	exp1 WT siCtrl	WT siCtrl	NT	2420	0.000000	1293.25	1.8712546	1.966433	0.9515984
WT	exp1	siPARP1	exp1 WT siPARP1	WT siPARP1	NT	2180	0.000000	824.50	2.6440267	2.608265	1.0137111
ALC1KO	exp1	siCtrl	exp1 ALC1KO siCtrl	ALC1KO siCtrl	NT	2428	0.000000	922.50	2.6319783	2.598642	1.0128285
ALC1KO	exp1	siPARP1	exp1 ALC1KO siPARP1	ALC1KO siPARP1	NT	2188	0.000000	848.00	2.5801887	2.618868	0.9852307
WT	exp2	siCtrl	exp2 WT siCtrl	WT siCtrl	NT	2072	0.000000	1049.00	1.9752145	1.966433	1.0044656
WT	exp2	siPARP1	exp2 WT siPARP1	WT siPARP1	NT	2504	0.000000	895.75	2.7954228	2.608265	1.0717558
ALC1KO	exp2	siCtrl	exp2 ALC1KO siCtrl	ALC1KO siCtrl	NT	3160	0.000000	1223.00	2.5838103	2.598642	0.9942926
ALC1KO	exp2	siPARP1	exp2 ALC1KO siPARP1	ALC1KO siPARP1	NT	3536	0.000000	1315.25	2.6884623	2.618868	1.0265743
WT	exp3	siCtrl	exp3 WT siCtrl	WT siCtrl	NT	1768	0.000000	861.25	2.0528302	1.966433	1.0439359
WT	exp3	siPARP1	exp3 WT siPARP1	WT siPARP1	NT	1888	0.000000	791.50	2.3853443	2.608265	0.9145331
ALC1KO	exp3	siCtrl	exp3 ALC1KO siCtrl	ALC1KO siCtrl	NT	2455	0.000000	951.50	2.5801366	2.598642	0.9928789
ALC1KO	exp3	siPARP1	exp3 ALC1KO siPARP1	ALC1KO siPARP1	NT	3222	0.000000	1245.00	2.5879518	2.618868	0.9881950
WT	exp1	siCtrl	exp1 WT siCtrl	WT siCtrl	H2O2_5uM	1528	0.698970	1293.25	1.1815194	1.966433	0.6008440
WT	exp1	siPARP1	exp1 WT siPARP1	WT siPARP1	H2O2_5uM	755	0.698970	824.50	0.9157065	2.608265	0.3510788
ALC1KO	exp1	siCtrl	exp1 ALC1KO siCtrl	ALC1KO siCtrl	H2O2_5uM	923	0.698970	922.50	1.0005420	2.598642	0.3850250
ALC1KO	exp1	siPARP1	exp1 ALC1KO siPARP1	ALC1KO siPARP1	H2O2_5uM	834	0.698970	848.00	0.9834906	2.618868	0.3755404
WT	exp2	siCtrl	exp2 WT siCtrl	WT siCtrl	H2O2_5uM	1332	0.698970	1049.00	1.2697807	1.966433	0.6457279
WT	exp2	siPARP1	exp2 WT siPARP1	WT siPARP1	H2O2_5uM	624	0.698970	895.75	0.6966229	2.608265	0.2670829
ALC1KO	exp2	siCtrl	exp2 ALC1KO siCtrl	ALC1KO siCtrl	H2O2_5uM	1222	0.698970	1223.00	0.9991823	2.598642	0.3845018
ALC1KO	exp2	siPARP1	exp2 ALC1KO siPARP1	ALC1KO siPARP1	H2O2_5uM	1296	0.698970	1315.25	0.9853640	2.618868	0.3762558
WT	exp3	siCtrl	exp3 WT siCtrl	WT siCtrl	H2O2_5uM	896	0.698970	861.25	1.0403483	1.966433	0.5290535
WT	exp3	siPARP1	exp3 WT siPARP1	WT siPARP1	H2O2_5uM	696	0.698970	791.50	0.8793430	2.608265	0.3371372
ALC1KO	exp3	siCtrl	exp3 ALC1KO siCtrl	ALC1KO siCtrl	H2O2_5uM	712	0.698970	951.50	0.7482922	2.598642	0.2879551
ALC1KO	exp3	siPARP1	exp3 ALC1KO siPARP1	ALC1KO siPARP1	H2O2_5uM	1224	0.698970	1245.00	0.9831325	2.618868	0.3754037
WT	exp1	siCtrl	exp1 WT siCtrl	WT siCtrl	H2O2_30uM	977	1.477121	1293.25	0.7554610	1.966433	0.3841784
WT	exp1	siPARP1	exp1 WT siPARP1	WT siPARP1	H2O2_30uM	303	1.477121	824.50	0.3674955	2.608265	0.1408965
ALC1KO	exp1	siCtrl	exp1 ALC1KO siCtrl	ALC1KO siCtrl	H2O2_30uM	244	1.477121	922.50	0.2644986	2.598642	0.1017834
ALC1KO	exp1	siPARP1	exp1 ALC1KO siPARP1	ALC1KO siPARP1	H2O2_30uM	296	1.477121	848.00	0.3490566	2.618868	0.1332853
WT	exp2	siCtrl	exp2 WT siCtrl	WT siCtrl	H2O2_30uM	598	1.477121	1049.00	0.5700667	1.966433	0.2898989
WT	exp2	siPARP1	exp2 WT siPARP1	WT siPARP1	H2O2_30uM	344	1.477121	895.75	0.3840357	2.608265	0.1472380
ALC1KO	exp2	siCtrl	exp2 ALC1KO siCtrl	ALC1KO siCtrl	H2O2_30uM	388	1.477121	1223.00	0.3172527	2.598642	0.1220840
ALC1KO	exp2	siPARP1	exp2 ALC1KO siPARP1	ALC1KO siPARP1	H2O2_30uM	396	1.477121	1315.25	0.3010834	2.618868	0.1149670
WT	exp3	siCtrl	exp3 WT siCtrl	WT siCtrl	H2O2_30uM	562	1.477121	861.25	0.6525399	1.966433	0.3318394
WT	exp3	siPARP1	exp3 WT siPARP1	WT siPARP1	H2O2_30uM	422	1.477121	791.50	0.5331649	2.608265	0.2044136
ALC1KO	exp3	siCtrl	exp3 ALC1KO siCtrl	ALC1KO siCtrl	H2O2_30uM	499	1.477121	951.50	0.5244351	2.598642	0.2018112
ALC1KO	exp3	siPARP1	exp3 ALC1KO siPARP1	ALC1KO siPARP1	H2O2_30uM	467	1.477121	1245.00	0.3751004	2.618868	0.1432300
WT	exp1	siCtrl	exp1 WT siCtrl	WT siCtrl	H2O2_50uM	248	1.698970	1293.25	0.1917649	1.966433	0.0975192
WT	exp1	siPARP1	exp1 WT siPARP1	WT siPARP1	H2O2_50uM	60	1.698970	824.50	0.0727714	2.608265	0.0279003
ALC1KO	exp1	siCtrl	exp1 ALC1KO siCtrl	ALC1KO siCtrl	H2O2_50uM	95	1.698970	922.50	0.1029810	2.598642	0.0396288
ALC1KO	exp1	siPARP1	exp1 ALC1KO siPARP1	ALC1KO siPARP1	H2O2_50uM	74	1.698970	848.00	0.0872642	2.618868	0.0333213
WT	exp2	siCtrl	exp2 WT siCtrl	WT siCtrl	H2O2_50uM	194	1.698970	1049.00	0.1849380	1.966433	0.0940475
WT	exp2	siPARP1	exp2 WT siPARP1	WT siPARP1	H2O2_50uM	111	1.698970	895.75	0.1239185	2.608265	0.0475099
ALC1KO	exp2	siCtrl	exp2 ALC1KO siCtrl	ALC1KO siCtrl	H2O2_50uM	122	1.698970	1223.00	0.0997547	2.598642	0.0383872
ALC1KO	exp2	siPARP1	exp2 ALC1KO siPARP1	ALC1KO siPARP1	H2O2_50uM	33	1.698970	1315.25	0.0250903	2.618868	0.0095806
WT	exp3	siCtrl	exp3 WT siCtrl	WT siCtrl	H2O2_50uM	219	1.698970	861.25	0.2542816	1.966433	0.1293111
WT	exp3	siPARP1	exp3 WT siPARP1	WT siPARP1	H2O2_50uM	160	1.698970	791.50	0.2021478	2.608265	0.0775028
ALC1KO	exp3	siCtrl	exp3 ALC1KO siCtrl	ALC1KO siCtrl	H2O2_50uM	140	1.698970	951.50	0.1471361	2.598642	0.0566204
ALC1KO	exp3	siPARP1	exp3 ALC1KO siPARP1	ALC1KO siPARP1	H2O2_50uM	67	1.698970	1245.00	0.0538153	2.618868	0.0205491

Plot Data

library(ggplot2)

# raw data
ggplot(dataset, aes(x=H2O2, y=Counts)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_smooth(method=lm, formula = y ~ poly(x,2), se=FALSE, aes(colour=siRNA)) +
        geom_point(aes(colour=siRNA, shape=Experiment), size=2) +        
        facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)") +
        scale_shape_manual(values=15:20) +
        scale_color_manual(values=c("#000000","#FF0000"))

# NormCounts Linear
ggplot(dataset, aes(x=H2O2, y=NormCounts, color=siRNA)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_point(aes(colour=siRNA), size=2) +        
        geom_smooth(method=lm, formula = y ~ x, se=FALSE) +
        facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)") +
        scale_color_manual(values=c("#000000","#FF0000"))

# NormCounts2 Linear
ggplot(dataset, aes(x=H2O2, y=NormCounts2, color=siRNA)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_point(aes(colour=siRNA), size=2) +        
        geom_smooth(method=lm, formula = y ~ x, se=FALSE) +
        facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)") +
        scale_color_manual(values=c("#000000","#FF0000"))

# NormCounts Quadratic
ggplot(dataset, aes(x=H2O2, y=NormCounts, color=siRNA)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_point(aes(colour=siRNA), size=2) +        
        geom_smooth(method=lm, formula = y ~ poly(x,2), se=FALSE) +
        facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)")+
        scale_color_manual(values=c("#000000","#FF0000"))

# NormCounts2 Quadratic
ggplot(dataset, aes(x=H2O2, y=NormCounts2, color=siRNA)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_point(aes(colour=siRNA), size=2) +        
        geom_smooth(method=lm, formula = y ~ poly(x,2), se=FALSE) +
        facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)") +
        scale_color_manual(values=c("#000000","#FF0000"))

# NormCounts Cubic
ggplot(dataset, aes(x=H2O2, y=NormCounts, color=siRNA)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_point(aes(colour=siRNA), size=2) +        
        geom_smooth(method=lm, formula = y ~ poly(x,3), se=FALSE) +
        facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)")+
        scale_color_manual(values=c("#000000","#FF0000"))

# NormCounts2 Cubic
ggplot(dataset, aes(x=H2O2, y=NormCounts2, color=siRNA)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_point(aes(colour=siRNA), size=2) +        
        geom_smooth(method=lm, formula = y ~ poly(x,3), se=FALSE) +
        facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)") +
        scale_color_manual(values=c("#000000","#FF0000"))

pdf("FigureS3_H2O2.pdf", width = 5, height = 4)

ggplot(dataset, aes(x=H2O2, y=NormCounts2)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_point(aes(colour = siRNA, shape = genotype), size=1.75) +
        geom_smooth(method=lm, formula = y ~ poly(x,2), se=TRUE, 
                    aes(group = GSID,colour = siRNA, linetype = genotype), fill='#DDDDDD', size=0.5) +
        #facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)") +
        ylab(label = "Normalized Counts") +
        scale_color_manual(values=c("#000000","#FF0000")) +
        guides(linetype = guide_legend(override.aes= list(color = "#555555"))) 

dev.off()

## quartz_off_screen 
##                 2

Models

library(MASS)
library(DHARMa)
library(lme4)
library(lmerTest)
library(bbmle)

Linear formula

fit1 <- lm(Counts ~ Experiment + H2O2*siRNA*genotype, data = dataset)
print(summary(fit1))

## 
## Call:
## lm(formula = Counts ~ Experiment + H2O2 * siRNA * genotype, data = dataset)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -697.48 -230.97    3.66  207.16  711.48 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                       2008.20     197.38  10.174 2.11e-12 ***
## Experimentexp2                     148.69     127.73   1.164   0.2516    
## Experimentexp3                      -9.75     127.73  -0.076   0.9396    
## H2O2                             -1018.50     155.34  -6.557 9.85e-08 ***
## siRNAsiPARP1                      -117.27     258.93  -0.453   0.6532    
## genotypeALC1KO                     361.79     258.93   1.397   0.1704    
## H2O2:siRNAsiPARP1                 -116.96     219.68  -0.532   0.5975    
## H2O2:genotypeALC1KO               -410.10     219.68  -1.867   0.0697 .  
## siRNAsiPARP1:genotypeALC1KO        423.11     366.19   1.155   0.2551    
## H2O2:siRNAsiPARP1:genotypeALC1KO   -91.64     310.68  -0.295   0.7696    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 361.3 on 38 degrees of freedom
## Multiple R-squared:  0.8871, Adjusted R-squared:  0.8603 
## F-statistic: 33.17 on 9 and 38 DF,  p-value: 2.56e-15

cat("AIC: ", AIC(fit1))

## AIC:  712.4085

simres <- simulateResiduals(fittedModel = fit1)
plot(simres)

fit2 <- lm(NormCounts ~ H2O2*siRNA*genotype, data = dataset)
print(summary(fit2))

## 
## Call:
## lm(formula = NormCounts ~ H2O2 * siRNA * genotype, data = dataset)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.66792 -0.08492  0.07217  0.16452  0.48643 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                        1.9297     0.1435  13.447  < 2e-16 ***
## H2O2                              -0.9597     0.1218  -7.882 1.14e-09 ***
## siRNAsiPARP1                       0.3793     0.2029   1.869   0.0690 .  
## genotypeALC1KO                     0.4093     0.2029   2.017   0.0505 .  
## H2O2:siRNAsiPARP1                 -0.3915     0.1722  -2.274   0.0284 *  
## H2O2:genotypeALC1KO               -0.4225     0.1722  -2.454   0.0186 *  
## siRNAsiPARP1:genotypeALC1KO       -0.3275     0.2870  -1.141   0.2606    
## H2O2:siRNAsiPARP1:genotypeALC1KO   0.3381     0.2435   1.388   0.1727    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2832 on 40 degrees of freedom
## Multiple R-squared:  0.9189, Adjusted R-squared:  0.9047 
## F-statistic: 64.75 on 7 and 40 DF,  p-value: < 2.2e-16

cat("AIC: ", AIC(fit2))

## AIC:  24.33749

simres <- simulateResiduals(fittedModel = fit2)
plot(simres)

fit3 <- lm(NormCounts2 ~ H2O2*siRNA*genotype, data = dataset)
print(summary(fit3))

## 
## Call:
## lm(formula = NormCounts2 ~ H2O2 * siRNA * genotype, data = dataset)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.25608 -0.04318  0.02865  0.07163  0.18649 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                       0.98133    0.05621  17.457  < 2e-16 ***
## H2O2                             -0.48804    0.04769 -10.233  9.9e-13 ***
## siRNAsiPARP1                     -0.09607    0.07950  -1.208    0.234    
## genotypeALC1KO                   -0.08123    0.07950  -1.022    0.313    
## H2O2:siRNAsiPARP1                -0.03001    0.06745  -0.445    0.659    
## H2O2:genotypeALC1KO              -0.04386    0.06745  -0.650    0.519    
## siRNAsiPARP1:genotypeALC1KO       0.10889    0.11243   0.969    0.339    
## H2O2:siRNAsiPARP1:genotypeALC1KO  0.01371    0.09539   0.144    0.886    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1109 on 40 degrees of freedom
## Multiple R-squared:  0.9246, Adjusted R-squared:  0.9114 
## F-statistic: 70.09 on 7 and 40 DF,  p-value: < 2.2e-16

cat("AIC: ", AIC(fit3))

## AIC:  -65.63346

simres <- simulateResiduals(fittedModel = fit3)
plot(simres)

fit4 <- lmer(Counts ~ H2O2*siRNA*genotype + (1|UID), data = dataset)
print(summary(fit4))

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Counts ~ H2O2 * siRNA * genotype + (1 | UID)
##    Data: dataset
## 
## REML criterion at convergence: 601.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0069 -0.4188  0.1341  0.4384  2.2311 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  UID      (Intercept)   3255    57.06  
##  Residual             127693   357.34  
## Number of obs: 48, groups:  UID, 12
## 
## Fixed effects:
##                                  Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)                       2054.52     184.07    35.29  11.161 3.97e-13
## H2O2                             -1018.50     153.65    32.00  -6.629 1.78e-07
## siRNAsiPARP1                      -117.27     260.32    35.29  -0.451   0.6551
## genotypeALC1KO                     361.79     260.32    35.29   1.390   0.1733
## H2O2:siRNAsiPARP1                 -116.96     217.29    32.00  -0.538   0.5941
## H2O2:genotypeALC1KO               -410.10     217.29    32.00  -1.887   0.0682
## siRNAsiPARP1:genotypeALC1KO        423.11     368.15    35.29   1.149   0.2582
## H2O2:siRNAsiPARP1:genotypeALC1KO   -91.64     307.30    32.00  -0.298   0.7675
##                                     
## (Intercept)                      ***
## H2O2                             ***
## siRNAsiPARP1                        
## genotypeALC1KO                      
## H2O2:siRNAsiPARP1                   
## H2O2:genotypeALC1KO              .  
## siRNAsiPARP1:genotypeALC1KO         
## H2O2:siRNAsiPARP1:genotypeALC1KO    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                (Intr) H2O2   sRNAsPARP1 gALC1K H2O2:sRNAPARP1 H2O2:A sRNAPARP1:
## H2O2           -0.809                                                          
## siRNAsPARP1    -0.707  0.572                                                   
## gntypALC1KO    -0.707  0.572  0.500                                            
## H2O2:sRNAPARP1  0.572 -0.707 -0.809     -0.404                                 
## H2O2:ALC1KO     0.572 -0.707 -0.404     -0.809  0.500                          
## sRNAPARP1:A     0.500 -0.404 -0.707     -0.707  0.572          0.572           
## H2O2:RNAPARP1: -0.404  0.500  0.572      0.572 -0.707         -0.707 -0.809

cat("AIC: ", AIC(fit4))

## AIC:  621.279

simres <- simulateResiduals(fittedModel = fit4)
plot(simres)

Quadratic formula

fit5 <- lm(Counts ~ Experiment + poly(H2O2,2)*siRNA*genotype, data = dataset)
print(summary(fit5))

## 
## Call:
## lm(formula = Counts ~ Experiment + poly(H2O2, 2) * siRNA * genotype, 
##     data = dataset)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -725.74 -158.61  -17.09  102.58  475.60 
## 
## Coefficients:
##                                            Estimate Std. Error t value Pr(>|t|)
## (Intercept)                                 1021.52      98.18  10.405 4.18e-12
## Experimentexp2                               148.69      98.18   1.514   0.1392
## Experimentexp3                                -9.75      98.18  -0.099   0.9215
## poly(H2O2, 2)1                             -4737.38     555.39  -8.530 5.80e-10
## poly(H2O2, 2)2                                61.12     555.39   0.110   0.9130
## siRNAsiPARP1                                -230.58     113.37  -2.034   0.0498
## genotypeALC1KO                               -35.50     113.37  -0.313   0.7561
## poly(H2O2, 2)1:siRNAsiPARP1                 -544.03     785.44  -0.693   0.4932
## poly(H2O2, 2)2:siRNAsiPARP1                 1647.83     785.44   2.098   0.0434
## poly(H2O2, 2)1:genotypeALC1KO              -1907.53     785.44  -2.429   0.0206
## poly(H2O2, 2)2:genotypeALC1KO               1754.59     785.44   2.234   0.0322
## siRNAsiPARP1:genotypeALC1KO                  334.33     160.33   2.085   0.0446
## poly(H2O2, 2)1:siRNAsiPARP1:genotypeALC1KO  -426.25    1110.78  -0.384   0.7036
## poly(H2O2, 2)2:siRNAsiPARP1:genotypeALC1KO -1694.41    1110.78  -1.525   0.1364
##                                               
## (Intercept)                                ***
## Experimentexp2                                
## Experimentexp3                                
## poly(H2O2, 2)1                             ***
## poly(H2O2, 2)2                                
## siRNAsiPARP1                               *  
## genotypeALC1KO                                
## poly(H2O2, 2)1:siRNAsiPARP1                   
## poly(H2O2, 2)2:siRNAsiPARP1                *  
## poly(H2O2, 2)1:genotypeALC1KO              *  
## poly(H2O2, 2)2:genotypeALC1KO              *  
## siRNAsiPARP1:genotypeALC1KO                *  
## poly(H2O2, 2)1:siRNAsiPARP1:genotypeALC1KO    
## poly(H2O2, 2)2:siRNAsiPARP1:genotypeALC1KO    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 277.7 on 34 degrees of freedom
## Multiple R-squared:  0.9403, Adjusted R-squared:  0.9175 
## F-statistic:  41.2 on 13 and 34 DF,  p-value: < 2.2e-16

cat("AIC: ", AIC(fit5))

## AIC:  689.8121

simres <- simulateResiduals(fittedModel = fit5)
plot(simres)

fit6 <- lm(NormCounts ~ poly(H2O2,2)*siRNA*genotype, data = dataset)
print(summary(fit6))

## 
## Call:
## lm(formula = NormCounts ~ poly(H2O2, 2) * siRNA * genotype, data = dataset)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.233442 -0.077535 -0.008586  0.088724  0.298066 
## 
## Coefficients:
##                                              Estimate Std. Error t value
## (Intercept)                                 1.000e+00  4.284e-02  23.343
## poly(H2O2, 2)1                             -4.464e+00  2.968e-01 -15.040
## poly(H2O2, 2)2                              1.093e-01  2.968e-01   0.368
## siRNAsiPARP1                                7.412e-17  6.058e-02   0.000
## genotypeALC1KO                             -5.228e-17  6.058e-02   0.000
## poly(H2O2, 2)1:siRNAsiPARP1                -1.821e+00  4.197e-01  -4.338
## poly(H2O2, 2)2:siRNAsiPARP1                 1.907e+00  4.197e-01   4.544
## poly(H2O2, 2)1:genotypeALC1KO              -1.965e+00  4.197e-01  -4.682
## poly(H2O2, 2)2:genotypeALC1KO               1.670e+00  4.197e-01   3.978
## siRNAsiPARP1:genotypeALC1KO                 9.821e-17  8.568e-02   0.000
## poly(H2O2, 2)1:siRNAsiPARP1:genotypeALC1KO  1.572e+00  5.936e-01   2.649
## poly(H2O2, 2)2:siRNAsiPARP1:genotypeALC1KO -2.133e+00  5.936e-01  -3.593
##                                            Pr(>|t|)    
## (Intercept)                                 < 2e-16 ***
## poly(H2O2, 2)1                              < 2e-16 ***
## poly(H2O2, 2)2                             0.714841    
## siRNAsiPARP1                               1.000000    
## genotypeALC1KO                             1.000000    
## poly(H2O2, 2)1:siRNAsiPARP1                0.000111 ***
## poly(H2O2, 2)2:siRNAsiPARP1                6.00e-05 ***
## poly(H2O2, 2)1:genotypeALC1KO              3.95e-05 ***
## poly(H2O2, 2)2:genotypeALC1KO              0.000322 ***
## siRNAsiPARP1:genotypeALC1KO                1.000000    
## poly(H2O2, 2)1:siRNAsiPARP1:genotypeALC1KO 0.011915 *  
## poly(H2O2, 2)2:siRNAsiPARP1:genotypeALC1KO 0.000970 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1484 on 36 degrees of freedom
## Multiple R-squared:   0.98,  Adjusted R-squared:  0.9738 
## F-statistic:   160 on 11 and 36 DF,  p-value: < 2.2e-16

cat("AIC: ", AIC(fit6))

## AIC:  -34.74272

simres <- simulateResiduals(fittedModel = fit6)
plot(simres)

fit7 <- lm(NormCounts2 ~ poly(H2O2,2)*siRNA*genotype, data = dataset)
print(summary(fit7))

## 
## Call:
## lm(formula = NormCounts2 ~ poly(H2O2, 2) * siRNA * genotype, 
##     data = dataset)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.098930 -0.035263 -0.003283  0.033948  0.125891 
## 
## Coefficients:
##                                            Estimate Std. Error t value Pr(>|t|)
## (Intercept)                                 0.50853    0.01779  28.581  < 2e-16
## poly(H2O2, 2)1                             -2.27004    0.12327 -18.415  < 2e-16
## poly(H2O2, 2)2                              0.05558    0.12327   0.451 0.654770
## siRNAsiPARP1                               -0.12514    0.02516  -4.973 1.63e-05
## genotypeALC1KO                             -0.12372    0.02516  -4.917 1.94e-05
## poly(H2O2, 2)1:siRNAsiPARP1                -0.13957    0.17433  -0.801 0.428600
## poly(H2O2, 2)2:siRNAsiPARP1                 0.71759    0.17433   4.116 0.000214
## poly(H2O2, 2)1:genotypeALC1KO              -0.20399    0.17433  -1.170 0.249633
## poly(H2O2, 2)2:genotypeALC1KO               0.62899    0.17433   3.608 0.000930
## siRNAsiPARP1:genotypeALC1KO                 0.12217    0.03559   3.433 0.001517
## poly(H2O2, 2)1:siRNAsiPARP1:genotypeALC1KO  0.06375    0.24654   0.259 0.797440
## poly(H2O2, 2)2:siRNAsiPARP1:genotypeALC1KO -0.80896    0.24654  -3.281 0.002301
##                                               
## (Intercept)                                ***
## poly(H2O2, 2)1                             ***
## poly(H2O2, 2)2                                
## siRNAsiPARP1                               ***
## genotypeALC1KO                             ***
## poly(H2O2, 2)1:siRNAsiPARP1                   
## poly(H2O2, 2)2:siRNAsiPARP1                ***
## poly(H2O2, 2)1:genotypeALC1KO                 
## poly(H2O2, 2)2:genotypeALC1KO              ***
## siRNAsiPARP1:genotypeALC1KO                ** 
## poly(H2O2, 2)1:siRNAsiPARP1:genotypeALC1KO    
## poly(H2O2, 2)2:siRNAsiPARP1:genotypeALC1KO ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06164 on 36 degrees of freedom
## Multiple R-squared:  0.9791, Adjusted R-squared:  0.9727 
## F-statistic:   153 on 11 and 36 DF,  p-value: < 2.2e-16

cat("AIC: ", AIC(fit7))

## AIC:  -119.0964

simres <- simulateResiduals(fittedModel = fit7)
plot(simres)

fit8 <- lmer(Counts ~ poly(H2O2,2)*siRNA*genotype + (1|UID), data = dataset)
print(summary(fit8))

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Counts ~ poly(H2O2, 2) * siRNA * genotype + (1 | UID)
##    Data: dataset
## 
## REML criterion at convergence: 505
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.44832 -0.62194  0.01954  0.46624  1.90599 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  UID      (Intercept) 19567    139.9   
##  Residual             62446    249.9   
## Number of obs: 48, groups:  UID, 12
## 
## Fixed effects:
##                                            Estimate Std. Error       df t value
## (Intercept)                                 1067.83     108.29     8.00   9.861
## poly(H2O2, 2)1                             -4737.38     499.78    28.00  -9.479
## poly(H2O2, 2)2                                61.12     499.78    28.00   0.122
## siRNAsiPARP1                                -230.58     153.14     8.00  -1.506
## genotypeALC1KO                               -35.50     153.14     8.00  -0.232
## poly(H2O2, 2)1:siRNAsiPARP1                 -544.03     706.80    28.00  -0.770
## poly(H2O2, 2)2:siRNAsiPARP1                 1647.83     706.80    28.00   2.331
## poly(H2O2, 2)1:genotypeALC1KO              -1907.53     706.80    28.00  -2.699
## poly(H2O2, 2)2:genotypeALC1KO               1754.59     706.80    28.00   2.482
## siRNAsiPARP1:genotypeALC1KO                  334.33     216.57     8.00   1.544
## poly(H2O2, 2)1:siRNAsiPARP1:genotypeALC1KO  -426.25     999.57    28.00  -0.426
## poly(H2O2, 2)2:siRNAsiPARP1:genotypeALC1KO -1694.41     999.57    28.00  -1.695
##                                            Pr(>|t|)    
## (Intercept)                                9.42e-06 ***
## poly(H2O2, 2)1                             3.09e-10 ***
## poly(H2O2, 2)2                               0.9035    
## siRNAsiPARP1                                 0.1706    
## genotypeALC1KO                               0.8225    
## poly(H2O2, 2)1:siRNAsiPARP1                  0.4479    
## poly(H2O2, 2)2:siRNAsiPARP1                  0.0272 *  
## poly(H2O2, 2)1:genotypeALC1KO                0.0117 *  
## poly(H2O2, 2)2:genotypeALC1KO                0.0193 *  
## siRNAsiPARP1:genotypeALC1KO                  0.1612    
## poly(H2O2, 2)1:siRNAsiPARP1:genotypeALC1KO   0.6731    
## poly(H2O2, 2)2:siRNAsiPARP1:genotypeALC1KO   0.1011    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                      (Intr) pl(H2O2,2)1 pl(H2O2,2)2 sRNAsPARP1 gALC1K
## pl(H2O2,2)1           0.000                                          
## pl(H2O2,2)2           0.000  0.000                                   
## siRNAsPARP1          -0.707  0.000       0.000                       
## gntypALC1KO          -0.707  0.000       0.000       0.500           
## pl(H2O2,2)1:RNAPARP1  0.000 -0.707       0.000       0.000      0.000
## pl(H2O2,2)2:RNAPARP1  0.000  0.000      -0.707       0.000      0.000
## p(H2O2,2)1:A          0.000 -0.707       0.000       0.000      0.000
## p(H2O2,2)2:A          0.000  0.000      -0.707       0.000      0.000
## sRNAPARP1:A           0.500  0.000       0.000      -0.707     -0.707
## p(H2O2,2)1:RNAPARP1:  0.000  0.500       0.000       0.000      0.000
## p(H2O2,2)2:RNAPARP1:  0.000  0.000       0.500       0.000      0.000
##                      pl(H2O2,2)1:RNAPARP1 pl(H2O2,2)2:RNAPARP1 p(H2O2,2)1:A
## pl(H2O2,2)1                                                                
## pl(H2O2,2)2                                                                
## siRNAsPARP1                                                                
## gntypALC1KO                                                                
## pl(H2O2,2)1:RNAPARP1                                                       
## pl(H2O2,2)2:RNAPARP1  0.000                                                
## p(H2O2,2)1:A          0.500                0.000                           
## p(H2O2,2)2:A          0.000                0.500                0.000      
## sRNAPARP1:A           0.000                0.000                0.000      
## p(H2O2,2)1:RNAPARP1: -0.707                0.000               -0.707      
## p(H2O2,2)2:RNAPARP1:  0.000               -0.707                0.000      
##                      p(H2O2,2)2:A sRNAPARP1: p(H2O2,2)1:RNAPARP1:
## pl(H2O2,2)1                                                      
## pl(H2O2,2)2                                                      
## siRNAsPARP1                                                      
## gntypALC1KO                                                      
## pl(H2O2,2)1:RNAPARP1                                             
## pl(H2O2,2)2:RNAPARP1                                             
## p(H2O2,2)1:A                                                     
## p(H2O2,2)2:A                                                     
## sRNAPARP1:A           0.000                                      
## p(H2O2,2)1:RNAPARP1:  0.000        0.000                         
## p(H2O2,2)2:RNAPARP1: -0.707        0.000      0.000

cat("AIC: ", AIC(fit8))

## AIC:  533.0262

simres <- simulateResiduals(fittedModel = fit8)
plot(simres)

Cubic formula

fit9 <- lm(Counts ~ Experiment + poly(H2O2,3)*siRNA*genotype, data = dataset)
print(summary(fit9))

## 
## Call:
## lm(formula = Counts ~ Experiment + poly(H2O2, 3) * siRNA * genotype, 
##     data = dataset)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -747.69 -119.62    6.48  112.67  451.62 
## 
## Coefficients:
##                                             Estimate Std. Error t value
## (Intercept)                                 1021.521     94.245  10.839
## Experimentexp2                               148.688     94.245   1.578
## Experimentexp3                                -9.750     94.245  -0.103
## poly(H2O2, 3)1                             -4737.379    533.132  -8.886
## poly(H2O2, 3)2                                61.116    533.132   0.115
## poly(H2O2, 3)3                              -743.243    533.132  -1.394
## siRNAsiPARP1                                -230.583    108.825  -2.119
## genotypeALC1KO                               -35.500    108.825  -0.326
## poly(H2O2, 3)1:siRNAsiPARP1                 -544.032    753.963  -0.722
## poly(H2O2, 3)2:siRNAsiPARP1                 1647.828    753.963   2.186
## poly(H2O2, 3)3:siRNAsiPARP1                    9.997    753.963   0.013
## poly(H2O2, 3)1:genotypeALC1KO              -1907.534    753.963  -2.530
## poly(H2O2, 3)2:genotypeALC1KO               1754.587    753.963   2.327
## poly(H2O2, 3)3:genotypeALC1KO                106.013    753.963   0.141
## siRNAsiPARP1:genotypeALC1KO                  334.333    153.902   2.172
## poly(H2O2, 3)1:siRNAsiPARP1:genotypeALC1KO  -426.247   1066.265  -0.400
## poly(H2O2, 3)2:siRNAsiPARP1:genotypeALC1KO -1694.412   1066.265  -1.589
## poly(H2O2, 3)3:siRNAsiPARP1:genotypeALC1KO   -54.356   1066.265  -0.051
##                                            Pr(>|t|)    
## (Intercept)                                6.78e-12 ***
## Experimentexp2                               0.1251    
## Experimentexp3                               0.9183    
## poly(H2O2, 3)1                             6.64e-10 ***
## poly(H2O2, 3)2                               0.9095    
## poly(H2O2, 3)3                               0.1735    
## siRNAsiPARP1                                 0.0425 *  
## genotypeALC1KO                               0.7465    
## poly(H2O2, 3)1:siRNAsiPARP1                  0.4761    
## poly(H2O2, 3)2:siRNAsiPARP1                  0.0368 *  
## poly(H2O2, 3)3:siRNAsiPARP1                  0.9895    
## poly(H2O2, 3)1:genotypeALC1KO                0.0169 *  
## poly(H2O2, 3)2:genotypeALC1KO                0.0269 *  
## poly(H2O2, 3)3:genotypeALC1KO                0.8891    
## siRNAsiPARP1:genotypeALC1KO                  0.0379 *  
## poly(H2O2, 3)1:siRNAsiPARP1:genotypeALC1KO   0.6922    
## poly(H2O2, 3)2:siRNAsiPARP1:genotypeALC1KO   0.1225    
## poly(H2O2, 3)3:siRNAsiPARP1:genotypeALC1KO   0.9597    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 266.6 on 30 degrees of freedom
## Multiple R-squared:  0.9515, Adjusted R-squared:  0.924 
## F-statistic:  34.6 on 17 and 30 DF,  p-value: 4.022e-15

cat("AIC: ", AIC(fit9))

## AIC:  687.8777

simres <- simulateResiduals(fittedModel = fit9)
plot(simres)

fit10 <- lm(NormCounts ~ poly(H2O2,3)*siRNA*genotype, data = dataset)
print(summary(fit10))

## 
## Call:
## lm(formula = NormCounts ~ poly(H2O2, 3) * siRNA * genotype, data = dataset)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.222920 -0.039175 -0.001219  0.045161  0.187158 
## 
## Coefficients:
##                                              Estimate Std. Error t value
## (Intercept)                                 1.000e+00  2.783e-02  35.936
## poly(H2O2, 3)1                             -4.464e+00  1.928e-01 -23.154
## poly(H2O2, 3)2                              1.093e-01  1.928e-01   0.567
## poly(H2O2, 3)3                             -6.838e-01  1.928e-01  -3.547
## siRNAsiPARP1                                1.666e-16  3.935e-02   0.000
## genotypeALC1KO                              1.302e-16  3.935e-02   0.000
## poly(H2O2, 3)1:siRNAsiPARP1                -1.821e+00  2.727e-01  -6.679
## poly(H2O2, 3)2:siRNAsiPARP1                 1.907e+00  2.727e-01   6.996
## poly(H2O2, 3)3:siRNAsiPARP1                -1.882e-01  2.727e-01  -0.690
## poly(H2O2, 3)1:genotypeALC1KO              -1.965e+00  2.727e-01  -7.208
## poly(H2O2, 3)2:genotypeALC1KO               1.670e+00  2.727e-01   6.124
## poly(H2O2, 3)3:genotypeALC1KO               5.081e-02  2.727e-01   0.186
## siRNAsiPARP1:genotypeALC1KO                -3.274e-16  5.565e-02   0.000
## poly(H2O2, 3)1:siRNAsiPARP1:genotypeALC1KO  1.572e+00  3.856e-01   4.078
## poly(H2O2, 3)2:siRNAsiPARP1:genotypeALC1KO -2.133e+00  3.856e-01  -5.531
## poly(H2O2, 3)3:siRNAsiPARP1:genotypeALC1KO  2.267e-01  3.856e-01   0.588
##                                            Pr(>|t|)    
## (Intercept)                                 < 2e-16 ***
## poly(H2O2, 3)1                              < 2e-16 ***
## poly(H2O2, 3)2                             0.574725    
## poly(H2O2, 3)3                             0.001227 ** 
## siRNAsiPARP1                               1.000000    
## genotypeALC1KO                             1.000000    
## poly(H2O2, 3)1:siRNAsiPARP1                1.54e-07 ***
## poly(H2O2, 3)2:siRNAsiPARP1                6.30e-08 ***
## poly(H2O2, 3)3:siRNAsiPARP1                0.495028    
## poly(H2O2, 3)1:genotypeALC1KO              3.47e-08 ***
## poly(H2O2, 3)2:genotypeALC1KO              7.59e-07 ***
## poly(H2O2, 3)3:genotypeALC1KO              0.853350    
## siRNAsiPARP1:genotypeALC1KO                1.000000    
## poly(H2O2, 3)1:siRNAsiPARP1:genotypeALC1KO 0.000281 ***
## poly(H2O2, 3)2:siRNAsiPARP1:genotypeALC1KO 4.23e-06 ***
## poly(H2O2, 3)3:siRNAsiPARP1:genotypeALC1KO 0.560695    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0964 on 32 degrees of freedom
## Multiple R-squared:  0.9925, Adjusted R-squared:  0.989 
## F-statistic: 281.6 on 15 and 32 DF,  p-value: < 2.2e-16

cat("AIC: ", AIC(fit10))

## AIC:  -73.81496

simres <- simulateResiduals(fittedModel = fit10)
plot(simres)

fit11 <- lm(NormCounts2 ~ poly(H2O2,3)*siRNA*genotype, data = dataset)
print(summary(fit11))

## 
## Call:
## lm(formula = NormCounts2 ~ poly(H2O2, 3) * siRNA * genotype, 
##     data = dataset)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.085467 -0.014959 -0.000465  0.019616  0.071756 
## 
## Coefficients:
##                                             Estimate Std. Error t value
## (Intercept)                                 0.508535   0.011511  44.179
## poly(H2O2, 3)1                             -2.270035   0.079750 -28.464
## poly(H2O2, 3)2                              0.055582   0.079750   0.697
## poly(H2O2, 3)3                             -0.347722   0.079750  -4.360
## siRNAsiPARP1                               -0.125138   0.016279  -7.687
## genotypeALC1KO                             -0.123719   0.016279  -7.600
## poly(H2O2, 3)1:siRNAsiPARP1                -0.139574   0.112783  -1.238
## poly(H2O2, 3)2:siRNAsiPARP1                 0.717592   0.112783   6.363
## poly(H2O2, 3)3:siRNAsiPARP1                 0.013414   0.112783   0.119
## poly(H2O2, 3)1:genotypeALC1KO              -0.203992   0.112783  -1.809
## poly(H2O2, 3)2:genotypeALC1KO               0.628990   0.112783   5.577
## poly(H2O2, 3)3:genotypeALC1KO               0.104147   0.112783   0.923
## siRNAsiPARP1:genotypeALC1KO                 0.122166   0.023022   5.307
## poly(H2O2, 3)1:siRNAsiPARP1:genotypeALC1KO  0.063748   0.159500   0.400
## poly(H2O2, 3)2:siRNAsiPARP1:genotypeALC1KO -0.808958   0.159500  -5.072
## poly(H2O2, 3)3:siRNAsiPARP1:genotypeALC1KO  0.003173   0.159500   0.020
##                                            Pr(>|t|)    
## (Intercept)                                 < 2e-16 ***
## poly(H2O2, 3)1                              < 2e-16 ***
## poly(H2O2, 3)2                             0.490862    
## poly(H2O2, 3)3                             0.000126 ***
## siRNAsiPARP1                               9.20e-09 ***
## genotypeALC1KO                             1.17e-08 ***
## poly(H2O2, 3)1:siRNAsiPARP1                0.224889    
## poly(H2O2, 3)2:siRNAsiPARP1                3.82e-07 ***
## poly(H2O2, 3)3:siRNAsiPARP1                0.906070    
## poly(H2O2, 3)1:genotypeALC1KO              0.079901 .  
## poly(H2O2, 3)2:genotypeALC1KO              3.71e-06 ***
## poly(H2O2, 3)3:genotypeALC1KO              0.362697    
## siRNAsiPARP1:genotypeALC1KO                8.14e-06 ***
## poly(H2O2, 3)1:siRNAsiPARP1:genotypeALC1KO 0.692050    
## poly(H2O2, 3)2:siRNAsiPARP1:genotypeALC1KO 1.61e-05 ***
## poly(H2O2, 3)3:siRNAsiPARP1:genotypeALC1KO 0.984250    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03987 on 32 degrees of freedom
## Multiple R-squared:  0.9922, Adjusted R-squared:  0.9886 
## F-statistic: 271.6 on 15 and 32 DF,  p-value: < 2.2e-16

cat("AIC: ", AIC(fit11))

## AIC:  -158.557

simres <- simulateResiduals(fittedModel = fit11)
plot(simres)

fit12 <- lmer(Counts ~ poly(H2O2,3)*siRNA*genotype + (1|UID), data = dataset)
print(summary(fit12))

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Counts ~ poly(H2O2, 3) * siRNA * genotype + (1 | UID)
##    Data: dataset
## 
## REML criterion at convergence: 439.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.67826 -0.45216 -0.01205  0.40378  1.92855 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  UID      (Intercept) 22071    148.6   
##  Residual             52429    229.0   
## Number of obs: 48, groups:  UID, 12
## 
## Fixed effects:
##                                             Estimate Std. Error        df
## (Intercept)                                 1067.833    108.288     8.000
## poly(H2O2, 3)1                             -4737.379    457.950    24.000
## poly(H2O2, 3)2                                61.116    457.950    24.000
## poly(H2O2, 3)3                              -743.243    457.950    24.000
## siRNAsiPARP1                                -230.583    153.142     8.000
## genotypeALC1KO                               -35.500    153.142     8.000
## poly(H2O2, 3)1:siRNAsiPARP1                 -544.032    647.639    24.000
## poly(H2O2, 3)2:siRNAsiPARP1                 1647.828    647.639    24.000
## poly(H2O2, 3)3:siRNAsiPARP1                    9.997    647.639    24.000
## poly(H2O2, 3)1:genotypeALC1KO              -1907.534    647.639    24.000
## poly(H2O2, 3)2:genotypeALC1KO               1754.587    647.639    24.000
## poly(H2O2, 3)3:genotypeALC1KO                106.013    647.639    24.000
## siRNAsiPARP1:genotypeALC1KO                  334.333    216.575     8.000
## poly(H2O2, 3)1:siRNAsiPARP1:genotypeALC1KO  -426.247    915.899    24.000
## poly(H2O2, 3)2:siRNAsiPARP1:genotypeALC1KO -1694.412    915.899    24.000
## poly(H2O2, 3)3:siRNAsiPARP1:genotypeALC1KO   -54.356    915.899    24.000
##                                            t value Pr(>|t|)    
## (Intercept)                                  9.861 9.42e-06 ***
## poly(H2O2, 3)1                             -10.345 2.53e-10 ***
## poly(H2O2, 3)2                               0.133  0.89495    
## poly(H2O2, 3)3                              -1.623  0.11766    
## siRNAsiPARP1                                -1.506  0.17057    
## genotypeALC1KO                              -0.232  0.82250    
## poly(H2O2, 3)1:siRNAsiPARP1                 -0.840  0.40919    
## poly(H2O2, 3)2:siRNAsiPARP1                  2.544  0.01780 *  
## poly(H2O2, 3)3:siRNAsiPARP1                  0.015  0.98781    
## poly(H2O2, 3)1:genotypeALC1KO               -2.945  0.00706 ** 
## poly(H2O2, 3)2:genotypeALC1KO                2.709  0.01225 *  
## poly(H2O2, 3)3:genotypeALC1KO                0.164  0.87134    
## siRNAsiPARP1:genotypeALC1KO                  1.544  0.16123    
## poly(H2O2, 3)1:siRNAsiPARP1:genotypeALC1KO  -0.465  0.64585    
## poly(H2O2, 3)2:siRNAsiPARP1:genotypeALC1KO  -1.850  0.07666 .  
## poly(H2O2, 3)3:siRNAsiPARP1:genotypeALC1KO  -0.059  0.95317    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("AIC: ", AIC(fit12))

## AIC:  475.7654

simres <- simulateResiduals(fittedModel = fit12)
plot(simres)

Compare Results

ICtab(fit1,fit2,fit3,fit4,
      fit5,fit6,fit7,fit8,
      fit9,fit10,fit11,fit12,
      base=T)

##       AIC    dAIC   df
## fit11 -158.6    0.0 17
## fit7  -119.1   39.5 13
## fit10  -73.8   84.7 17
## fit3   -65.6   92.9 9 
## fit6   -34.7  123.8 13
## fit2    24.3  182.9 9 
## fit12  475.8  634.3 18
## fit8   533.0  691.6 14
## fit4   621.3  779.8 10
## fit9   687.9  846.4 19
## fit5   689.8  848.4 15
## fit1   712.4  871.0 11

Final Result

fit <- fit11

output <- coef(summary(fit))
output <- output[grep("H2O2", rownames(output)),]


rownames(output) <- gsub("poly\\(|, [1-3]\\)","", rownames(output) )
rownames(output) <- gsub("genotype",  paste0(" ",levels(dataset$genotype)[1], " vs. "), rownames(output))
rownames(output)[!(grepl("vs", rownames(output)))] <- paste(rownames(output)[!(grepl("vs", rownames(output)))], levels(dataset$genotype)[1],  sep = " in " )

rownames(output) <- gsub("siRNA",  paste0(" ",levels(dataset$siRNA)[1], " vs. "), rownames(output))
rownames(output)[!(grepl("vs.*vs| in ", rownames(output)))] <- paste(rownames(output)[!(grepl("vs.*vs| in ", rownames(output)))], levels(dataset$siRNA)[1],  sep = " in " )

rownames(output)[!(grepl("vs", rownames(output)))] <- paste(rownames(output)[!(grepl("vs", rownames(output)))], levels(dataset$siRNA)[1],  sep = " " )


# suggested result table
kable(output, row.names = T)

	Estimate	Std. Error	t value	Pr(>|t|)
H2O21 in WT siCtrl	-2.2700352	0.0797498	-28.4644590	0.0000000
H2O22 in WT siCtrl	0.0555822	0.0797498	0.6969573	0.4908621
H2O23 in WT siCtrl	-0.3477220	0.0797498	-4.3601612	0.0001261
H2O21: siCtrl vs. siPARP1 in WT	-0.1395743	0.1127833	-1.2375448	0.2248890
H2O22: siCtrl vs. siPARP1 in WT	0.7175920	0.1127833	6.3625752	0.0000004
H2O23: siCtrl vs. siPARP1 in WT	0.0134139	0.1127833	0.1189354	0.9060700
H2O21: WT vs. ALC1KO in siCtrl	-0.2039924	0.1127833	-1.8087113	0.0799010
H2O22: WT vs. ALC1KO in siCtrl	0.6289899	0.1127833	5.5769793	0.0000037
H2O23: WT vs. ALC1KO in siCtrl	0.1041468	0.1127833	0.9234246	0.3626969
H2O21: siCtrl vs. siPARP1: WT vs. ALC1KO	0.0637482	0.1594996	0.3996763	0.6920503
H2O22: siCtrl vs. siPARP1: WT vs. ALC1KO	-0.8089581	0.1594996	-5.0718495	0.0000161
H2O23: siCtrl vs. siPARP1: WT vs. ALC1KO	0.0031734	0.1594996	0.0198959	0.9842500

write.table(output, file = "FigureS3_H2O2_Stats_Ref_WT.txt", quote = F, sep = "\t", row.names = T, col.names = NA)

# re-fit with ALC1KO reference
dataset$genotype <- relevel(dataset$genotype, ref = "ALC1KO")


fit <- lm(NormCounts2 ~ poly(H2O2,3)*siRNA*genotype, data = dataset)

output <- coef(summary(fit))
output <- output[grep("H2O2", rownames(output)),]


rownames(output) <- gsub("poly\\(|, [1-3]\\)","", rownames(output) )
rownames(output) <- gsub("genotype",  paste0(" ",levels(dataset$genotype)[1], " vs. "), rownames(output))
rownames(output)[!(grepl("vs", rownames(output)))] <- paste(rownames(output)[!(grepl("vs", rownames(output)))], levels(dataset$genotype)[1],  sep = " in " )

rownames(output) <- gsub("siRNA",  paste0(" ",levels(dataset$siRNA)[1], " vs. "), rownames(output))
rownames(output)[!(grepl("vs.*vs| in ", rownames(output)))] <- paste(rownames(output)[!(grepl("vs.*vs| in ", rownames(output)))], levels(dataset$siRNA)[1],  sep = " in " )

rownames(output)[!(grepl("vs", rownames(output)))] <- paste(rownames(output)[!(grepl("vs", rownames(output)))], levels(dataset$siRNA)[1],  sep = " " )


# suggested result table
kable(output, row.names = T)

	Estimate	Std. Error	t value	Pr(>|t|)
H2O21 in ALC1KO siCtrl	-2.4740275	0.0797498	-31.0223630	0.0000000
H2O22 in ALC1KO siCtrl	0.6845721	0.0797498	8.5839971	0.0000000
H2O23 in ALC1KO siCtrl	-0.2435752	0.0797498	-3.0542416	0.0045206
H2O21: siCtrl vs. siPARP1 in ALC1KO	-0.0758261	0.1127833	-0.6723172	0.5062076
H2O22: siCtrl vs. siPARP1 in ALC1KO	-0.0913661	0.1127833	-0.8101032	0.4238645
H2O23: siCtrl vs. siPARP1 in ALC1KO	0.0165873	0.1127833	0.1470725	0.8839975
H2O21: ALC1KO vs. WT in siCtrl	0.2039924	0.1127833	1.8087113	0.0799010
H2O22: ALC1KO vs. WT in siCtrl	-0.6289899	0.1127833	-5.5769793	0.0000037
H2O23: ALC1KO vs. WT in siCtrl	-0.1041468	0.1127833	-0.9234246	0.3626969
H2O21: siCtrl vs. siPARP1: ALC1KO vs. WT	-0.0637482	0.1594996	-0.3996763	0.6920503
H2O22: siCtrl vs. siPARP1: ALC1KO vs. WT	0.8089581	0.1594996	5.0718495	0.0000161
H2O23: siCtrl vs. siPARP1: ALC1KO vs. WT	-0.0031734	0.1594996	-0.0198959	0.9842500

write.table(output, file = "FigureS3_H2O2_Stats_Ref_ALC1KO.txt", quote = F, sep = "\t", row.names = T, col.names = NA)

Anova

fit11a <- lm(NormCounts2 ~ poly(H2O2,3)*siRNA*genotype, data = dataset)
fit11b <- lm(NormCounts2 ~ poly(H2O2,3)*siRNA+genotype, data = dataset)

# anova table
anova(fit11a, fit11b)

## Analysis of Variance Table
## 
## Model 1: NormCounts2 ~ poly(H2O2, 3) * siRNA * genotype
## Model 2: NormCounts2 ~ poly(H2O2, 3) * siRNA + genotype
##   Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
## 1     32 0.05088                                  
## 2     39 0.15961 -7  -0.10873 9.7692 1.861e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

fit11c <- lm(NormCounts2 ~ poly(H2O2,3)*genotype*siRNA, data = dataset)
fit11d <- lm(NormCounts2 ~ poly(H2O2,3)*genotype+siRNA, data = dataset)

# anova table
anova(fit11c, fit11d)

## Analysis of Variance Table
## 
## Model 1: NormCounts2 ~ poly(H2O2, 3) * genotype * siRNA
## Model 2: NormCounts2 ~ poly(H2O2, 3) * genotype + siRNA
##   Res.Df     RSS Df Sum of Sq      F   Pr(>F)    
## 1     32 0.05088                                 
## 2     39 0.16428 -7  -0.11339 10.188 1.21e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Links

• https://bioconductor.org/packages/release/bioc/vignettes/CFAssay/inst/doc/cfassay.pdf
• http://www.flutterbys.com.au/stats/course.html
• https://bookdown.org/roback/bookdown-bysh/
• https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html
• https://en.wikipedia.org/wiki/Akaike_information_criterion